# 井壁稳定性解释
import math

from seaman.core import JsonUtil
from wellstability.stability_common import StabilityParams




class ShaftResult:
    sigma_xx: float
    sigma_yy: float
    sigma_zz: float
    tau_xy: float
    tau_yz: float
    tau_xz: float


# 井壁围岩应力
def get_shaft_lining_result(params: StabilityParams) -> ShaftResult:
    result = ShaftResult()
    max_horizontal_ground_stress: float = params.max_horizontal_ground_stress_gradient * (params.depth / 1000)
    min_horizontal_ground_stress: float = params.min_horizontal_ground_stress_gradient * (params.depth / 1000)
    vertical_ground_stress: float = params.vertical_ground_stress_gradient * (params.depth / 1000)
    # 求解sigma_xx
    sigma_xx = (max_horizontal_ground_stress * math.pow(math.cos(math.radians(params.stability_azimuth)), 2)
                + min_horizontal_ground_stress * math.pow(math.sin(math.radians(params.stability_azimuth)), 2)) * math.pow(math.cos(math.radians(params.stability_inclination)), 2) + vertical_ground_stress * math.pow(math.sin(math.radians(params.stability_inclination)), 2)
    result.sigma_xx = sigma_xx
    # 求解sigma_yy
    sigma_yy = max_horizontal_ground_stress * math.pow(math.sin(math.radians(params.stability_azimuth)), 2) + min_horizontal_ground_stress * math.pow(math.cos(math.radians(params.stability_azimuth)), 2)
    result.sigma_yy = sigma_yy
    # 求解sigma_zz
    sigma_zz = (max_horizontal_ground_stress * math.pow(math.cos(math.radians(params.stability_azimuth)), 2)
                + min_horizontal_ground_stress * math.pow(math.sin(math.radians(params.stability_azimuth)), 2)) * math.pow(math.sin(math.radians(params.stability_inclination)), 2) + vertical_ground_stress * math.pow(math.cos(math.radians(params.stability_inclination)), 2)
    result.sigma_zz = sigma_zz
    # 求解tau_xy
    tau_xy = (min_horizontal_ground_stress - max_horizontal_ground_stress) * 0.5 * math.sin(math.radians(2 * params.stability_azimuth)) * math.cos(math.radians(params.stability_inclination))
    result.tau_xy = tau_xy
    # 求解tau_yz
    tau_yz = (min_horizontal_ground_stress - max_horizontal_ground_stress) * 0.5 * math.sin(math.radians(2 * params.stability_azimuth)) * math.sin(math.radians(params.stability_inclination))
    result.tau_yz = tau_yz
    # 求解tau_xz
    tau_xz = (max_horizontal_ground_stress * math.pow(math.cos(math.radians(params.stability_azimuth)), 2)
              + min_horizontal_ground_stress * math.pow(math.sin(math.radians(params.stability_azimuth)), 2)
              - vertical_ground_stress) * 0.5 * math.sin(math.radians(2 * params.stability_inclination))
    result.tau_xz = tau_xz
    return result


class RoundResult:
    optimal_theta: float
    z: float


# 井周角计算-------从0---180迭代计算 z取最小时 拿到theta值
def get_round_angle_result(params: StabilityParams) -> RoundResult:
    result = RoundResult()
    temp_one: ShaftResult = get_shaft_lining_result(params)
    sigma_xx = temp_one.sigma_xx
    sigma_yy = temp_one.sigma_yy
    tau_xy = temp_one.tau_xy
    max_theta = 180
    step_size = 0.1  # 调整步长为 0.1 度
    min_z = float('inf')  # 用于跟踪最小的 z 值
    optimal_theta = 0  # 用于跟踪对应最小 z 值的 theta
    theta = 0
    while theta <= max_theta:
        z = 2 * (sigma_xx - sigma_yy) * math.cos(2 * math.radians(theta)) + 4 * tau_xy * math.sin(
            2 * math.radians(theta))
        if z < min_z:
            min_z = z
            optimal_theta = theta
        theta += step_size
    result.z = min_z
    result.optimal_theta = optimal_theta
    return result


class MiddleResult:
    R: float
    K: float


# 失稳准则： D.P准则
def get_middle_result(params: StabilityParams) -> MiddleResult:
    result = MiddleResult()
    factor_a: float = math.pow(math.sin(math.radians(params.internal_friction_angle)), 2)
    R = (math.sqrt(3) * math.sin(math.radians(params.internal_friction_angle))) / (3 * math.sqrt(3 + factor_a))
    result.R = R
    K = (math.sqrt(3) * params.cohesion * math.cos(math.radians(params.internal_friction_angle))) / math.sqrt(3 + factor_a)
    result.K = K

    return result

# 计算Pm临界液体压力  临界密度
class CriticalResult:
    sigma_r: float
    sigma_theta: float
    sigma_z: float
    tau_theta_z: float
    tau_rz: float
    tau_r_theta: float
    sigma_one: float
    sigma_two: float
    sigma_three: float
    J_two: float
    I_one: float
    p_m: float


#  计算Pm  临界密度
def get_critical_result(params: StabilityParams) -> CriticalResult:
    result = CriticalResult()
    in_situ_pore_pressure: float = params.original_liquid_density * 9.8 * (params.depth / 1000)
    temp_one: ShaftResult = get_shaft_lining_result(params)
    sigma_xx = temp_one.sigma_xx
    sigma_yy = temp_one.sigma_yy
    sigma_zz = temp_one.sigma_zz
    tau_xy = temp_one.tau_xy
    tau_yz = temp_one.tau_yz
    tau_xz = temp_one.tau_xz
    temp_two = get_round_angle_result(params)
    theta = temp_two.optimal_theta
    temp_three = get_middle_result(params)
    R = temp_three.R
    K = temp_three.K
    p_m: float = 1.0
    # 设置迭代的要求
    step_size = 0.1  # 迭代步长为 0.1
    # 定义停止条件
    tolerance = 1e-6  # 容差值
    max_iterations = 1000  # 最大迭代次数
    iteration = 0
    # theta = 86.3
    # 求解三个带有p_m的公式
    while True:

        # *******************************
        # 计算sigma_r
        sigma_r = p_m - params.osmotic_coefficient * params.void_ratio * (p_m - in_situ_pore_pressure)
        # 计算sigma_theta
        factor_a: float = (params.formation_stress_factor * (1 - 2 * params.poisson_ratio)) / (1 - params.poisson_ratio)
        sigma_theta = sigma_xx + sigma_yy - 2 * (sigma_xx - sigma_yy) * math.cos(math.radians(2 * theta)) - 4 * tau_xy * math.sin(math.radians(2 * theta)) - p_m + params.osmotic_coefficient * (p_m - in_situ_pore_pressure) * (factor_a - params.void_ratio)
        # 计算sigma_z
        sigma_z = sigma_zz - (2 * params.poisson_ratio * (sigma_xx - sigma_yy) * math.cos(math.radians(2 * theta)) - 4 * params.poisson_ratio * tau_xy * math.sin(math.radians(2 * theta))) + params.osmotic_coefficient * (p_m - in_situ_pore_pressure) * (factor_a * 0.5 - params.void_ratio)
        # *******************************

        # 计算tau_theta_z
        tau_theta_z = 2 * tau_yz * math.cos(math.radians(theta)) - 2 * tau_xz * math.sin(math.radians(theta))
        # 计算tau_rz, tau_r_theta
        tau_rz = 0
        tau_r_theta = 0
        # 计算sigma_one
        sigma_one = sigma_r - params.formation_stress_factor * in_situ_pore_pressure
        # 计算sigma_two
        factor_b: float = (sigma_theta + sigma_z) / 2
        factor_c: float = math.pow((sigma_theta - sigma_z), 2) + 4 *math.pow(tau_theta_z, 2)
        factor_d: float = params.formation_stress_factor * in_situ_pore_pressure
        sigma_two = (factor_b + math.sqrt(factor_c) / 2) - factor_d
        # 计算sigma_three
        sigma_three = (factor_b - math.sqrt(factor_c) / 2) - factor_d
        # 计算J_two
        factor_e: float = math.pow((sigma_one - sigma_two), 2)
        factor_f: float = math.pow((sigma_two - sigma_three), 2)
        factor_g: float = math.pow((sigma_three - sigma_one), 2)
        J_two = (factor_e + factor_f + factor_g) / 6
        # 计算I_one
        I_one = sigma_one + sigma_two + sigma_three
        y_result = math.sqrt(J_two) - R * I_one - K

        # 更新 p_m 的值
        p_m += step_size
        # 检查停止条件
        if y_result < tolerance:

            break

        iteration += 1
    result.p_m = p_m
    result.sigma_r = sigma_r
    result.sigma_theta = sigma_theta
    result.sigma_z = sigma_z
    result.tau_theta_z = tau_theta_z
    result.tau_rz = tau_rz
    result.tau_r_theta = tau_r_theta
    result.sigma_one = sigma_one
    result.sigma_two = sigma_two
    result.sigma_three = sigma_three
    result.J_two = J_two
    result.I_one = I_one

    return result


class DrillingFluidDensityResult:
    rh0: float
    xb: float


# 临界钻井液密度
def get_drilling_fluid_density_result(params: StabilityParams) -> DrillingFluidDensityResult:
    result = DrillingFluidDensityResult()
    temp_one = get_critical_result(params)
    p_m = temp_one.p_m
    # p_m: float = 27.4
    rh0 = p_m / params.depth * 100
    result.rh0 = rh0
    xb: float = params.drilling_fluid_density / rh0
    result.xb = xb
    return result


class StrainResult:
    N: float


def get_strain_result(params: StabilityParams) -> StrainResult:
    result = StrainResult()
    max_horizontal_ground_stress: float = params.max_horizontal_ground_stress_gradient * (params.depth / 1000)
    min_horizontal_ground_stress: float = params.min_horizontal_ground_stress_gradient * (params.depth / 1000)
    vertical_ground_stress: float = params.vertical_ground_stress_gradient * (params.depth / 1000)
    N = ((3 * 2 * params.borehole_radius) / params.yang_modulus_elasticity) * ((max_horizontal_ground_stress + min_horizontal_ground_stress) + 2 * (max_horizontal_ground_stress - min_horizontal_ground_stress))
    result.N = N
    return result


class RunResult:
    rh0: float
    xb: float
    N: float


def run(params: StabilityParams) -> RunResult:
    result = RunResult()
    temp_one = get_drilling_fluid_density_result(params)
    rh0 = temp_one.rh0
    result.rh0 = rh0
    xb = temp_one.xb
    result.xb = xb
    temp_two = get_strain_result(params)
    N = temp_two.N
    result.N = N
    return result


if __name__ == '__main__':
    params = StabilityParams()
    params.depth = 1000
    params.cohesion = 2
    params.internal_friction_angle = 34
    params.formation_stress_factor = 0.7
    params.poisson_ratio = 0.2
    params.max_horizontal_ground_stress_gradient = 18.3
    params.min_horizontal_ground_stress_gradient = 17.8
    params.vertical_ground_stress_gradient = 20.5
    params.original_liquid_density = 1
    params.borehole_radius = 10.8
    params.void_ratio = 0.02
    params.osmotic_coefficient = 1
    params.yang_modulus_elasticity = 20000
    params.drilling_fluid_density = 1.52
    params.stability_inclination = 31.433
    params.stability_azimuth = 98.499

    # params.depth = 1000
    # params.cohesion = 2
    # params.internal_friction_angle = 34
    # params.formation_stress_factor = 0.7
    # params.poisson_ratio = 0.2
    # params.max_horizontal_ground_stress = 18.8
    # params.min_horizontal_ground_stress = 17
    # params.vertical_ground_stress = 20.5
    # params.in_situ_pore_pressure = 9.9
    # params.borehole_radius = 10.8
    # params.void_ratio = 0.02
    # params.osmotic_coefficient = 1
    # params.yang_modulus_elasticity = 20000
    # params.drilling_fluid_density = 1.46
    # params.stability_inclination = 30
    # params.stability_azimuth = 46
    shaft_result = get_shaft_lining_result(params)
    round_result = get_round_angle_result(params)
    middle_result = get_middle_result(params)
    Four_result = get_critical_result(params)
    Six_result = get_strain_result(params)
    Five_result = get_drilling_fluid_density_result(params)
    run_ruselt = run(params)

    print('===当前压力值为===')
    print(JsonUtil.obj2str(shaft_result))
    print('===当前井周角为===')
    print(JsonUtil.obj2str(round_result))
    print('===当前失稳准则为===')
    print(JsonUtil.obj2str(middle_result))
    print('===当前临界液体压力和临界密度为===')
    print(JsonUtil.obj2str(Four_result))
    print("===当前临界钻井液密度为===")
    print(JsonUtil.obj2str(Five_result))
    print('===当前计算最大变形为===')
    print(JsonUtil.obj2str(Six_result))
    print('===最后结果总览===')
    print(JsonUtil.obj2str(run_ruselt))


